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We discuss thermalization of isolated quantum systems by using the Husimi-Wehrl entropy eval¬ 
uated in the semiclassical treatment. The Husimi-Wehrl entropy is the Wehrl entropy obtained 
by using the Husimi function for the phase space distribution. The time evolution of the Husimi 
function is given by smearing the Wigner function, whose time evolution is obtained in the semi¬ 
classical approximation. We show the efficiency and usefulness of this semiclassical treatment 
in describing entropy production of a couple of quantum mechanical systems, whose classi¬ 
cal counter systems are known to be chaotic. We propose two methods to evaluate the time 
evolution of the Husimi-Wehrl entropy, the test-particle method and the two-step Monte-Carlo 
method. We demonstrate the characteristics of the two methods by numerical calculations, and 
show that the simultaneous application of the two methods ensures the reliability of the results 
of the Husimi-Wehrl entropy at a given time. 
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1. Introduction 

Thermalization process or entropy creation of isolated quantum systems is a long-standing issue, but 
not well understood problem. Relevant systems include the early universe where the transition from a 
vacuum state to a thermalized state occurs at the end of cosmic inflation, and the QCD matter created 
in the initial stage of relativistic heavy-ion collisions where thermal matter should be formed in a 
rather short time. It is known that both systems are well described in semiclassical approximation, 
and moreover a chaotic behavior of the classical limit may play some role in the entropy production. 
The present paper is concerned with the entropy production of an isolated quantum system for which 
the semiclassical approximation is valid and the classical counter part may show a chaotic behavior. 

To describe entropy in a pure quantum system, one may of course adopt the von Neumann 
entropy 111] as quantum mechanical entropy given by 


^vN = - Tr [p log p] , 


( 1 ) 


where p is the density matrix. For a pure state, however, p is idempotent, p^ = p, implying that the 
eigen value of p is 0 or 1, and the von Neumann entropy is zero. Even if we start from a mixed 
state, the time evolution described by a unitary operator will never lead to entropy growth. On the 
other hand, the entropy production in a rarefied gas composed of classical or quantum mechanical 
particles can be well described by an analog of the H function of Boltzmann given in terms of the 
distribution function/(q,p): 


S = - 


d^qd^p 


f{q,p) log f{q,p). 


( 2 ) 


It is noteworthy that a phase-space description is desirable for making classical-quantum correspon¬ 
dence clear, and even natural when the semiclassical approximation is valid. The standard method 
for such a description is to use the celebrated Wigner function 0], which is defined as a Wigner 
transform of the density matrix: The Wigner function fw{Q,p) can be regarded as a quasi phase- 
space distribution function. The use of the Wigner function as the phase space distribution function 
= fw{<l^p) in Eq. ([Hi), however, has essential drawbacks: Eirst, the Wigner function is, 
actually, not a genuine distribution function; f„ can be negative, which prevents us to calculate the 
entropy density according to Eq. Q. Second, the entropy defined by Eq. Q given in terms of the 
Wigner function fw does not grow in time, because the Wigner transform only gives an equivalent 
description of the quantum system in terms of, say, the q- or p-representation ylldLoilalTlj. Some 
coarse graining of the phase space is needed to describe an entropy production. 

In a classical chaotic system, two adjacent points in the phase space depart from each other expo¬ 
nentially in time. If available phase space volume is limited, the exponentially diffusing classical 
trajectories have to be folded in a complicated manner in the phase space. After a certain time start¬ 
ing from a localized phase space cell, a given phase space cell {27rh)^ consists of the mixture of 
trajectories stemming from the initially occupied localized cell and vacant regions not yet visited. 
Since we cannot distinguish the phase space points in a cell due to the uncertainty principle, it is 
reasonable to define a phase space distribution as a smeared or coarse-grained function over the 
phase space cell. 

We adopt the Husimi function /«(<?, p) Isl] as such a coarse-grained distribution function, which is 
defined as the expectation value of the density matrix with respect to a coherent state |^;). It is readily 
shown that /«(<?, p) is semi-positive definite, /„ > 0, and a coarse-grained function of the Wigner 
function, as will be shown in a later section. It is shown that the Husimi function faithfully 
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describes the characteristic properties of the underlying classical system, and has been utilized to 


identify the chaotic remnants in quantum systems [Sl llO, lUl ll2l]. Thus a natural candidate of the 


quantum mechanical entropy is given by Q with f (q,p) being substituted by the Husimi function 
fH{Q,p)- This entropy was introduced by Wehrl 11311 and may be called the Wehrl entropy, although 
he himself called it the classical entropy and failed in identifyin g th e distribution function fniQ^p) 
with the Husimi function: Such an identification was made later llldll . We refer to the Wehrl entropy 
obtained by using the Husimi function as the Husimi-Wehrl (HW) entropy jlSll . 


•Shw = — 


d^qd 


D 


P 


(27rh) 


D 


fniq.P) log fniq^p) 


( 3 ) 


It is worth mentioning that the HW entropy can be a good measure for a quantum entanglement 


of a system including quantum optical systems lIlhL Il7h . For a one-dimensional case, there is a 
minimum of 5hw = 1 ilsUlQIl . in contrast to the von Neumann entropy, which takes 5vn = 0 in 
the ground state. It is also shown that the HW entropy takes a value close to the von Neumann 
entropy at high temperature, and its growth rate coincides with the Kolmogorov-Sinai entropy for 
the one-dimensional inverted harmonic oscillator lll5n . 

A direct evaluation of the HW entropy for a quantum system is a kind of challenge even for the 
system with a few degrees of freedom because it involves a large-dimensional integral over the phase 
space even apart from the cumbersome calculation of the logarithm with precision. Nevertheless 
the HW entropy and its time evolution have been calculated for some quantum systems i20L 12111 . 
The equation of motion (EOM) of the Husimi function is given in ilOll . which contains a term of 
the order h, and thus has a more complicated form than that of the Wigner function even in the 
semiclassical approximation; see below. To solve the complicated EOM of the Husimi function, a 
test-particle method was proposed by Tsai and Muller 12111 . where the evolution of the test particles 
are determined to reproduce some of the moments. 

As already mentioned, the semiclassical approximation is suitable to reveal the effect of the 
chaotic nature of the classical counter part. It is noteworthy that the time evolution of the Wigner 
function in the semiclassical approximation where the 0{h'^) terms are ignored is readily obtained 
by solving the classical Hamilton equation; quantum mechanical information such as the uncertainty 
relation is encoded in the initial Wigner function, provided that it is given as the Wigner transform 
of the quantum density matrix. 

The time evolution of the Husimi function is given by smearing the time-evolved Wigner func¬ 
tion obtained in the semiclassical approximation. This is the method we adopt in this article. We 
shall show its efficiency and usefulness in describing entropy production using a couple of quantum 
mechanical systems whose respective classical counter systems are known to be chaotic. We pro¬ 
pose two methods to evaluate the time evolution of the Husimi-Wehrl entropy. One is an adaptation 
of the usual test-particle method without recourse to the moments of the distribution function. The 
other is a sequential application of Monte-Carlo integration, which we call the two-step Monte-Carlo 
method. We shall demonstrate the characteristics of the two methods by numerical calculations, and 
show that the simultaneous application of the two methods ensures the reliability of the results of the 
HW-entropy’s time evolution. It should be noted that these two methods are, in principle, applicable 
to systems with large degrees of freedom such as quantum field fheories. 

The paper is organized as follows. In Sec. ID we summarize some basic ingredienfs of fhe Wigner 
and Husimi functions fogefher wifh fhe HW enfropy. In Sec. [3l we infroduce fhe fwo numerical 
mefhods fo evaluafe fhe HW enfropy in an efficienf way. In Sec. HI fhe quanfum mechanical models 
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are introduced and numerical results of the Husimi-Wehrl entropy are shown. The final section is 
devoted to a brief summary and concluding remarks. 

2. Wigner function, Husimi function, and Husimi-Wehrl entropy 

In this secti( 
and Husimi 

While we introduce Wigner and Husimi functions in one-dimensional quantum mechanics in 
Subsec. |2H]and|2]2j extension to multi-dimensional cases is straightforward. 


in, we briefly review quantum mechanical phase space distribution functions, Wigner O 
functions, and the phase space expression of the entropy, Husimi-Wehrl entropy ilSll . 


2.1. Wigner and Husimi functions 

The Wigner function 10] is defined as a Wigner transform of the density matrix 
fw{q,p,t) =Pw{q,p,t) 

= J dpe-^P^/\q+^\pit)\q-^). (4) 

While the Wigner function fw{q,p) can be regarded as a quasi phase space distribution function and 
provides intuitive picture of the phase space dynamics, it is not semi-positive definite and hence we 
cannot regard fw{q,p) as the phase space probability density. 

In order to overcome the above drawbacks of the Wigner function, Husimi introduced a Gaussian 
smeared Wigner function 1^, known as the Husimi function, 

fH{q,p) = I ^ ( 5 ) 

where A is an arbitrary width parameter that gives the smearing manner in the phase space. 

The Husimi function is defined also as the expectation value of the density matrix with respect to 
a coherent state | z): 

fH{q,p) = {z\p\z), z = {Aq + ip)/s/2hA, ( 6 ) 


for a one-dimensional case with A being an arbitrary constant. Here the coherent state is given by 

1^) = a = {Aq + ip)/V2hK, (7) 


where |0) is the ground state; a|0) = 0. It is readily shown that fniq-iP) is semi-positive definite, 
/h > 0 by using Eq. dh]); = \{z\f)\^ > 0 for a pure state \ f), and > 0 for 

a mixed state specified by the density matrix p = Yhi Wi\'4>i) > 0). 

The Husimi function fniq^p) serves as the probability density to observe the phase space variables 
{q, p) under a minimum wave packet \z), and is now semi-positive definite, /„ > 0. Compared with 
the Wigner function, the Husimi function is smooth and the peak of the Husimi function often 


appears around the expectation value of the position and momentum HlU. 


2.2. Time evolution in semiclassical approximation 

The equation of motion (EOM) for the Wigner function /jy is obtained from the Wigner transform 
of the von Neumann equation for the density matrix, dp/dt = [H, p]/ih. By applying the Wigner 
transform of the operator product, {AB)^ = exp{ih{^p — Vp^q)/2)H,y lESB], com¬ 
mutators are replaced by Poisson brackets as [A, B]^/ih = {A, i?}pB -h 0{h?). Thus the EOM for 
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fw is given in terms of the Wigner transform of the Hamiltonian H as 

dU 


dt 


= {H^,U}p^ + 0{h^) 


( 8 ) 


The Wigner transform Hy^ of a Hamiltonian with the form of if = p^/2m + U{q) does not change 
its form. We note that the 0{h?) term in dUl is proportional to the third derivative of Hy^ or U. Thus 
the EOM dUl without the 0{h^) term turns out to he exact for some simple models such as a (an 
inverted) harmonic oscillator. 

The semiclassical EOM for fyy is given hy retaining the terms up to 0{h) in Eq. dUl, which reads 

(9fiw (9/, 


dfyy dHyy Ofyy 


= 0 


(9) 


dt dp dq dq dp 

We remark that the semiclassical EOM is exact for the linear systems mentioned above. Equation 
da asserts that fyy is constant along the classical trajectory: Eet us see this. Eet {q{t; q), p{t;p)) is a 
solution of the classical EOM, i.e., Hamilton’s equation; 

dq _dHyy dp _ dHyy 

dt dp ' dt dq ^ 

with an initial condition (g(0) = q,p{0) = p). Then we have for fyy{q{t] q), p{t;p), t), 


( 10 ) 


DU 

Dt 


dfyy dqdfyy dpdfy 

dt ^ ^ ^ 


= 0 


( 11 ) 


dt dq dt dp 

which implies that is time-independent; fvy{q{t]q), p{t',p), t) = fyy{q, p, 0). Accordingly we 
have 


fw{Q, P, t) = fw{q{-t-,q), p{-t;p), 0). 


( 12 ) 


Thus we can obtain the semiclassical time evolution of the Wigner function by solving the classi¬ 
cal equation of motion. Note that the quantum mechanical effects are taken into account through 
the distribution of the initial value in the phase space encoded in the Wigner function /»(<?, p, 0) 
constructed from the initial density matrix. 

It is worth mentioning that the exact analytical solution of the time evolution of for some linear 
systems including a (stable) harmonic oscillator potential Is, 15|, an inverted (unstable) harmonic 


rUi 


oscillator potential lIlSL 12311 and an external potential 111511 can be obtained. Then even the analytic 
form of the Husimi function /^(q, p, t) for these systems are readily obtained ilSll by the Gaussian 
smearing of fw{Q,Pi which is easy to perform analytically. 

We note here that one may obtain the time evolution of the Husimi function /^(q, p, t) by solving 
the EOM for fniq, P, t), which involves terms proportional to h, and thus has a more complicated 
structure than that for fyy (g, p, t) even in the semiclassical approximation ilOll . If one sticks to solve 
the EOM for directly, some numerical method would be necessary. A test-particle method is 
adopted as such a numerical method by Tsai and Muller 12IH . where the time evolution of test parti¬ 
cles are determined so as to reproduce some of moments. We remark that there are some ambiguities 
in such an approach inherent in the moment method. 

In this work, we do not adopt this direct method for obtaining the time evolution of the Husimi 
function /^(g, p, t). We take advantage of the fact that the EOM of the Wigner function fyy{q, p, t) 
in this regime is given simply by solving the classical EOM, and obtain /^(q, p,t) by the Gaus¬ 
sian smearing of thus obtained fyy{q, p, t). This strategy should be workable and natural when the 
semiclassical approximation is meaningful. The remaining task that we have to do for obtaining the 
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Husimi function is just the multi-dimensional integrations over the phase space with the Gaussian 
kernel for the smearing, which should he feasible hy standard methods such as the Monte-Carlo 
integration. 


2.3. Husimi-Wehrl entropy 

Since the Wigner function f„ is merely the Weyl transform of the density matrix, any ohservahle 
is calculable in terms of fw in principle, and it is also the case with the Husimi function /„. A 
drawback of the is that it can have negative values, and hence is not suitable for the calculation 
of entropy. As is mentioned in Introduction and the previous subsection, the Husimi function is, 
in contrast, a semi-positive definite coarse-grained phase space distribution function smeared by 
a minimum wave packet, and hence a good candidate for the phase space distribution f{q,p) to 
evaluate the entropy of a quantum system, as the H function of Boltzmann in the classical system, 
Eq. (|2]), or equivalently the Husimi-Wehrl entropy given in Eq. (l3]l lll3n . 

An explicit form of the HW entropy in terms of the Wigner function is given by substituting the 
D-dimensional extension of Eq. (I5]l into Eq. (I3]l, 


X log 


pD i/jD ir 
a q a p -a o- 


(irh) 


D 


^-Aig-g'r/n-ip-p'r/An t) 


(13) 


One may now recognize some difficulty of the numerical evaluation of the HW entropy: It involves 
repeated numerical integrations over the multi-dimensional phase space, and in particular one of 
them appears as an argument of logarithm, which turns out to be quite problematic in the Monte- 
Carlo integration. 


3. Numerical methods to analyze the semiclassical time evolution of Husimi-Wehrl 
entropy 

Here, two numerical methods are introduced to calculate the time dependence of the HW entropy as 
given by the Gaussian smearing of the Wigner function obtained in the semiclassical approximation. 
Both methods are based on an adaptation of the Monte-Carlo integration over the phase-space. We 
call the two methods the test-particle (TP) and two-step Monte-Carlo (tsMC) methods, respectively. 
In this section, we deal with the D-dimensional system described by the Hamiltonian H = H{q, p), 
where q and p denote the D-dimensional vector, respectively, i.e., q = (gi, q 2 , ■ ■ ■ ,9d) and p = 
(pi, P2,... ,Pd)- 


3.1. Test-particle method 


In the test-particle method 1241251126112711 . the Wigner function is represented as a sum of the delta 
functions. 


19 N-r-p 

fw{q, p, t) =—^ — qi{t))3^{p-Pi{t )), 


A^tp ^ 


( 14 ) 


with the initial function 


2=1 


fry fy\D ^TP 

fw{q,p,o) = qi{^))3^{p-Pi{^)), 




i=l 
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where A^tp is the total number of the test particles, and their coordinates are given hy {qi{t), Pi{t)). 
The initial distrihution of the test particles ((?j(0), Pi{0)) {i = 1, 2,..., D) is chosen so as to 
well sample that of fw{Q, P, 0): Hence A^tp is called the sampling number. The time evolution 
of the coordinates {qi{t),pi{t)) is determined by the EOM for fw{q, p, t), which is reduced to the 
canonical equation of motion, 

dqi _ dH^ dpi _ dH^ ("151 

in the semiclassical approximation. 

For the test-particle representation of the Wigner function Eq. (fT4l) . the Husimi function is readily 
expressed as 


Ntp 


fH{q,P,t) = 


Nj'p 


E- 

2=1 


-Aiq-q.{t)r/h-ip-p,it)r/Ah 


(16) 


It is noteworthy that the Husimi function here is a smooth function in contrast to the corresponding 
Wigner function in Eq. (fT4l) . 

Inserting the Wigner function (fT4l) into Eq. (fTH) . the HW entropy in the test-particle method is 
given as. 


c.(TP) 


r d^qd^p A[q-q,[t)Y/n-(p-p,(t)Y/Ah i f „ f) 
Ntp^J {7rh)D ® ^og U[q,P,t) 


(17) 


Now note that the integral over {q,p)i for each i has a support only around the positions of the 
test particles {qi{t), Pi{t)) due to the Gaussian function, and then we can effectively perform the 
Monte-Carlo integration as follows; By generating a set of random numbers {Q, P)i with standard 
deviations of ^Jh/2K and \/%Kj2, Monte-Carlo sampling point {q,p)i for each i is obtained as 
{q^p)i = (Q) P)i + {di^Pi)- Thus we reach the formula to be used in the actual evaluation of the 
HW entropy in the test-particle method: 


o(TP) 

‘^HW 


_ 


1 

Wmc^tp 


A^mc A^tp 


I] I] log 


r,D 

j=i 


(18) 


where the amount of the sample number of (Q, P)i is denoted by A^mc- 


3.2. Two-step Monte-Carlo method 

The second method is a direct Monte-Carlo evaluation of the multi-dimensional integrals. We rewrite 
Eq. (fTH) as 


„(tsMC) 


r d^Qd^P AQpn-pyAH f d^qd^p 

J {ph)D J {27rh)D 


X log 


d^Q'd^P' 


g-A(Q')Vft-(P')VAft 


fwiq + Q + Q' tP + P + P'P) 


• out ’ in 

— “ log fwidk + Qk + Q’ljPk P Pk + -P/j t) 

out T H IIT T H 

k=l L 1=1 

= — (log {fw{q + Q + Q',p + P + P' 


QPqp 


(19) 
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where {Q^, P^) and (Q;, -P/) are Gaussian random numbers for the Monte-Carlo (MC) integration 
to compute Husimi function fniQ,?)- For the (g,p)-integration, we generate MC samples {q',p') 
at f = 0 according to the initial distribution, and obtain the corresponding phase space sample 
points {q{q',p', t),p{q', p', t)) at t by solving the canonical equation of motion. Under the semi- 
classical approximation, is constant and the Jacobian is unity along the classical trajectory, 
J {q{t), p{t)/q' {0) , p' (0)) = 1. Then we can replace the integral over (q,p) in the first line of Eq. (fT^ 
with the integral at f = 0 by using the initial distribution and the Liouville theorem as, 

f d^qd^p ^ ^ 

J {27rh)D fw{q,P,t)g{q,p) 

= / fwiq',p',0)g{q{q',p',t),p{q',p',t)) , (20) 

where {q',p') are the phase space coordinates at f = 0, and {q{q',p',t),p{q',p',t)) are those at t 
evolved from {q',p'). 

The Wigner function at t in the log in Eq. (fT^ can be obtained by the trace back of the trajectory 
from f to f = 0 as shown in Eq. dT^ . Equation (fT^ contains an MC integral of a function obtained 
by an MC integral; we first generate {q' ,p') at f = 0 according to the distribution fwiq^P, 0) and 
(Q, P) as Gaussian random numbers, and then perform the MC integral in the log by generating 
MC samples (Q', P'). We call this procedure two-step Monte-Carlo (tsMC). 

In the following sections, we show the characteristic properties of the two methods and 
demonstrate numerically how they work using two-dimensional quantum-mechanical systems. 


4. Numerical calculation of Husimi-Wehrl entropy in quantum Yang-Mills model 


In this section, we show the numerical results of the HW entropy in “quantum Yang-Mills system” 
1 2811 ■ obtained by the two distinct methods, TP and tsMC methods. 


4.1. Model Hamiltonian and setup of initial condition 
The Hamiltonian of the system is given by 

+P 2 ) + (21) 

We have restricted ourselves to the two-dimensional case here. The name, “quantum Yang-Mills 
(qYM)”, is originated from the fact that the spatially uniform Yang-Mills system is reduced to a 
(0 -I- 1)-dimensional system, i.e., a quantum mechanical system, and its Hamiltonian is just given 
by Eq. (EU). 

We adopt the initial condition given by a minimal wave packet centered at [qi,q 2 ,Pi,P 2 ) = 

( 0 , 0 , 10 , 10 ), 


fwiPi,P2,qi,q2,t = 0) = 


( 22 ) 


This initial condition is also adopted in Ref. 12111 . 

In the following, we show numerical results calculated by using the TP and tsMC methods. We 
show the results in the unit with m = 1 and h= 1, and take A = 1 for the wave packet width. In 
the case of A / 1, the smearing Gaussian is not symmetric in p and q directions. But the results do 
not change qualitatively. We have confirmed that the results with A = 0.1 and 10 are qualitatively 
the same as those with A = 1. 







Fig. 1 Time dependence of the HW entropy by using TP method in qYM, with A^tp = 
100,1000, 5000 and 15000, and A^mc = 500. The arrow shows how the calculated HW entropy 
changes as A^tp increases. 


4.2. Numerical results with TP method 

First, we show the numerical results of the HW entropy in the qYM system calculated in the TP 
method using Eq. (fT^ . 

Figure[T]shows the time evolution of the HW entropy calculated in the TP method with the follow¬ 
ing test-particle numbers, A^tp = 100,1000,5000 and 15000. The MC sample number is taken to 
be A^mc = 500. The statistical errors are estimated for A^mc samples from a standard deviation. We 
note that the calculated HW entropy at each t tends to increase along with increasing A^TP^ which 
is an artifact due to the small number of the test particles A^tp and discussed later. Apart from tiny 
fluctuations, all the calculation show that the HW entropy first increases in time with a small oscil¬ 
latory behavior being accompanied; its local maxima are seen around t ~ 0.5 and 1.7. We note that 
a similar behavior is also seen in Ref. Kill . 

Entropy evaluated by the TP method has a (unphysical) maximum depending on A^tp> which 
causes apparent saturation at large t in Fig. [T] In fact, when the system is chaotic and the phase 
space volume is very large, all the test particles will be so separated from each other in the phase 
space at later time that only the i = j terms in Eq. (fT^ will remain. In this limiting case, the HW 
entropy as given in (fTSl) is evaluated as follows; 

\ 


4 


(TP) 

HW 


log 


^ A^xp / 

= D — D log 2 -|- log Atp 


-NQyh-pyNh 


QP 


(23) 


which gives the inevitable upper limit of . In Appendix A, we examine the HW entropy of an 
inverted harmonic oscillator, for which 5 hw can be calculated analytically and is found to increase 
permanently. At later times, S'hw is underestimated with small Ntp values because of the upper 
limit discussed above. By comparison, S'hw at early times is calculated precisely in the TP method, 
as long as Atp is large enough for 5 hw to converge. 

From the above argument, S'Hw(f) would be obtained reliably as an extrapolated value in the 
limit of A^xp —)• oo. The extrapolation should be made in the A^xp range, where the limiting value 
is larger than the HW entropy to be obtained. The limiting values are ~ 
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Fig. 2 Time dependence of HW entropy calculated by using tsMC method. 


for A^tp = 100,1000,5000 and 15000, respectively. The large-t values found in Fig. [Dare close to 
these limiting values for smaller A^tp> i-e-> IVrp = 100 and 1000. Thus we see that the saturation 
behavior seen for smaller values of A^tp may be an artifact of the TP method. In contrast, the large- 
t values for A^tp = 5000 and 15000 in Fig. [Dare well below the limiting values (9.1 and 10.2), 
found free from the above mentioned artifact, and can be used to obtain the extrapolated value at 
A^tp —>• as discussed later in Subsec. lT4l Thus we conclude that the entropy production of the 
“quantum Yang-Mills” system can be well described with the use of HW entropy as calculated with 
the TP method with sufficiently large number of the test particles. 

4.3. Numerical results with tsMC method 

Next, we show the numerical results of the HW entropy in qYM in the tsMC method using the 
formula Eq. (fT^ . 

Figure |2] shows the time evolution of the HW entropy calculated in the tsMC method with the 
sample numbers Ni^ = 1200, 2400,4800 and 12000. Nout is taken to be the same as A^in. The errors 
attached to Shw in the present figure is estimated only for the Monte-Carlo integrals outside of log 
in Eq. (fT^ . and those from the integral inside the log is not taken into account, which causes an 
additional systematic error. 

We see that the larger the value of A^in, the smaller the HW entropy, which is an opposite depen¬ 
dence on the sample number to that in the TP method. Nevertheless the gross behavior in the time 
evolution of the HW entropy is quite similar in the two methods apart from the tiny fluctuations; 
After showing an oscillatory behavior in a first short period, it increases in a monotonous way 
and its growth rate decreases gradually. More quantitative comparison of the two methods will be 
presented in the next subsection. 

4.4. Comparison of the two methods 

Figure |3] shows the HW entropy at f = 10 as a function of Atp (-^in) in the TP (tsMC) method. 
We fit a linear function f{t) = at + b to the calculated 5Hw(i) data in the range 10 — Ai < f < 
10 -|- Af (Ai = 1), and adopt f{t= 10) as the HW entropy value at i = 10. This procedure provides 
a smoother curve and reduces the errors coming from fluctuations compared to directly using the 
raw data. 
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0 2 4 6 8 10 12 14 

Nxp, Ni„ (/1000) 


Fig. 3 HW entropy in qYM at t = 10 as a function of A^tp (-^Mc)^ and its extrapolation to 
infinitely large A^tp (-^mc) in the TP (tsMC) method. Filled circles (squares) show TP (tsMC) 
results, and the solid (dashed) line shows a fit function to TP (tsMC) results. The dotted line is 
limiting value given hy Eq. (|2^ . The shaded areas show the extrapolated value in the limit of 

A^tp, A^mc —^ oo. 


The HW entropy in the TP method becomes larger with increasing N^p as already mentioned; At 
t = 10, 5hw — 5.1 for A^xp = 100 and 5 hw — 8-7 for Ntp = 15000. We also show the fit results 
to the data for larger samples, say Ntp > 5000, with a fit function, 

= ( 2 '^) 

The extrapolated value to A^tp — oo is 9.19 ±0.10. When we use other fit functions such as 
f{N) = a — b/{N/c ± 1) and f{N) = a — b/N ± c/N^, the fit results have differences with a 
standard deviation of 0.16, which should he considered as a systematic error. Thus the HW entropy 
in the TP method is obtained as 

= 10) = 9.19 ± 0.10 (stat.) ± 0.16 (syst.) . (25) 

With increasing A^in, the HW entropy calculated in the tsMC method decreases, which is an oppo¬ 
site behavior to that in the TP method as noted before. At f = 10, S'hw — 13.2 for A^in = 1200 and 
Shw — 9.5 for A^in = 12000. We also show the fit results to the data. We adopt Eq. (1^ for the fit 
function. Erom the fit results, the HW entropy in the tsMC method is found to be 

= 10) = 9.01 ± 0.21 (stat.) ± 0.06 (syst.) , (26) 

where the central value and the statistical error are obtained from the fit using Eq. d24l) . and the 
systematic error is evaluated from the fits using several fit functions as done in the TP method. 

4.5. Discussions 

The time evolution of the HW entropies obtained in the TP and tsMC methods shows a similar 
behavior with each other: The HW entropy increases with an oscillatory behavior in the early stage, 
then shows a monotonous increase with a decreasing rate. The HW entropy at each t in the TP 
method increases along with ATP^ while it decreases with increasing in the tsMC method. 
Thus we can guess that the real value of the HW entropy lies between the results in the TP and 
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Fig. 4 Time dependence of HW entropy by using TP method in modified qYM. 


tsMC methods. Actually, the extrapolated values att = 10, = 10) = 9.19 ± 0.10 ± 0.16 at 

A^tp —>• oo and = 10) = 9.01 ± 0.21 ± 0.06 at Ajn —)• oo in the TP and tsMC methods 

respectively, are consistent with each other within the error. These results are also in agreement with 
that in Ref. 121]. 

These two methods, TP and tsMC methods, give consistent results after N ^ oo extrapolation. 
On the other hand, with finite number of A^tp and A^in, they could give seemingly inconsistent 
results depending on the dynamics. We here have a deeper look at this issue. In the tsMC method, 
the entropy seems to keep increasing even for the later time, in contrast to the results in the TP 
method with finite IVtp and in Ref. 112IH . The discrepancy may come from the special shape of the 
potential: there are two flat directions in the potential for the qYM system, although the width of 
them tends to shrink at large distances. Then, the classical trajectory can keep growing along the 
flat direction, which would cause an unlimited spreading of the Husimi function and a permanent 
increase of the HW entropy calculated in the semiclassical approximation. (In the case of the TP 
methods, there exists limiting value of the HW entropy depending on N'tp, which gives rise to the 
apparent saturation of S at large t. ) By comparison, it is shown that the exact energy spectra of 
the qYM are all discrete ones, because of the shrinking width leading to an increase of the kinetic 
energy due to the uncertainty relation, although the volume of {{p, q)\H{p, q) < E} is infinite 12911 . 
Note that the discrete spectra implies that the wave functions of the energy eigen states are all bound. 
Thus the corresponding Husimi function would not have a support at the infinite distance due to the 
quantum effect, and the HW entropy may not show the ever increasing behavior but have a saturated 
value. This plausible conjecture can only be confirmed by a full quanfum calculafion beyond fhe 
semiclassical approximation. Such a calculafion is beyond fhe scope of fhe presenf work and will be 
leff as a fufure work. Insfead, we shall fake anofher model, which is a modified version of fhe qYM 
one free from flat directions in its potential. 


5. Modified quantum Yang-Mills model 

Let us consider the model in which quartic potential terms are added to the qYM Hamiltonian; 


P\ 


P\ 




( 27 ) 
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Fig. 5 Time dependence of HW entropy by using tsMC method in modified qYM. 



Fig. 6 HW entropy in mqYM at f = 10 as a function of A^tp (-^mc)> and its extrapolation to 
infinitely large Wtp (-^mc) in the TP (tsMC) method. Filled circles (squares) show TP (tsMC) 
results, and the solid (dashed) line shows a fit function to TP (tsMC) results. The dotted line is 
limiting value given by Eq. (|2^ . The shaded areas show the extrapolated value in the limit of 

WtP, ^ OO. 


We call the system “modified quantum Yang-Mills (mqYM)”. The system is studied in Ref. il lLll2ll 
with < 0 in the context of chaos. It is apparent that there is no flat direction in the potential due 
to the quartic terms. We take = I and e = 0.1 in the Hamiltonian, Eq. (ITTl) . The mqYM system 
is found to be integrable with ej = 1,1/3 and oo llZBOll . Our choice oiejg^ = 0.1 is well apart 
from the integrable region. Since e is not very large, the HW entropy shows a similar behavior to 
that in qYM at early times, as shown later. 

In this section, we shall calculate the HW entropy of the mqYM system in the TP and tsMC 
methods . The analyses are carried out in a similar way to those for the qYM system. 

In Figs, m and [51 we show the time evolution of the HW entropy in mqYM calculated using the 
TP (Mtp = 500,1000, 5000 and 15000 with Nmc = 500) and tsMC (A^i„ = 600,1200,2400 and 
12000) methods, respectively. Nout is taken to be the same as Win for tsMC. 


13 















The distribution function in Eq. (l2^ is used as the initial condition, and the statistical errors are 
estimated for A^mc (Mn) samples from a standard deviation in the TP (tsMC) method, as in the 
qYM cases. 

Both of the calculated results show that the HW entropy first increases with an oscillatory behavior 
and tends to saturate at later times, t > 6. The later-time 5 hw values depend on the sample number, 
A^tp and A^in; With increasing Ntp (A^in)^ the HW entropy increases (decreases) in the TP (tsMC) 
method. These are the features also found in qYM. By comparison, it should be noted that there 
seems to be saturation of 5 hw both in the TP and tsMC methods in mqYM, in contrast to qYM. 
This may be originated from the finite phase space volume where the Husimi function has a support. 

In Fig.O we show the HW entropy at f = 10 as a function of A^tp or A^in. We fit a linear function 
to calculated S'Hw(f) results in the range 9 < f < 11, and adopt f{t = 10) as the HW entropy value 
at f = 10. In the TP method, = 10) ~ 6.4 and 7.5 for Ntp = 500 and 15000, respectively. 

In tsMC, we find = 10) ~ 9.4 and 7.7 for Win = 600 and 12000, respectively. 

The extrapolated values of 5 hw at Ntp —)• oo and Win ^ oo are found to be 

= 10) = 7.61 ± 0.01(stat.) ± 0.03(syst.) , (28) 

= 10) = 7.53 ± 0.01(stat.) ± 0.04(syst.) , (29) 

in the TP and tsMC methods, respectively. The central values and the statistical errors are obtained 
from the fit using Eq. m, and the systematic error is evaluated from the fits using several fit 
functions. These two values are consistent with each other within the error. 

The observation shows that the two methods, tsMC and TP, are especially effective for such a 
potential which bounds Husimi function in finite region. Thus, we are confident of the validity of 
the two methods in the mqYM system. 


6. Summary 

We have discussed entropy creation in isolated quantum systems by using the Husimi-Wehrl entropy 
evaluated in a semiclassical treatment. The semiclassical treatment is known to be useful in some of 
the systems such as the inflation in early universe and the early stage of relativistic heavy ion colli¬ 
sions. These systems are expected to bear instabilities and/or chaoticities in their classical counter 
systems, then the smearing of the phase space distribution by the minimal wave packet causes the 
entropy production in terms of the Wehrl entropy or the H function of Boltzmann even in isolated 
quantum systems. This is nothing but the Husimi-Wehrl entropy, the Wehrl entropy obtained by 
using the Gaussian smeared Wigner function (Husimi function) for the phase space distribution. 

The semiclassical time evolution of the Husimi function is given by solving a classical equation 
of motion and smearing with a Gaussian packet. Combining this semiclassical treatment with the 
Monte-Carlo numerical integral technique, we have developed two methods, the test-particle (TP) 
method and the two-step Monte Carlo (tsMC) method. We have applied these two methods to 
quantum mechanical systems in two dimensions, the quantum Yang-Mills (qYM) and the modi¬ 
fied quantum Yang-Mills (mqYM) systems. The classical counter systems of these are known to be 
chaotic. We have demonstrated that the Husimi-Wehrl entropy obtained in the TP (tsMC) method 
approaches the converged value from below (from above) with an increasing sample number, then 
we can guess the true value of HW entropy. We have further found that the results of the TP and tsMC 
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methods in the infinite sampling number limit are consistent within the error. Therefore, the simul¬ 
taneous application of the two methods ensures the reliahility of the results of the Husimi-Wehrl 
entropy at a given time. 

The extension of our methods to a multidimensional system is straightforward. We expect that 
these methods are useful in systems with many degrees of freedom such as the quantum field fheory. 
These mefhods are, in principle, applicable fo higher-dimensional problems, and we have confirmed 
fhaf fhey acfually work in fhree and four dimensional systems. In higher dimensions, we need much 
more Monfe-Carlo samples fo obfain sfafisfically reliable resulfs, and if would be necessary fo make 
some approximafions for practical purposes. Work in fhis direction is in progress. 
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A. HW entropy in inverted harmonics 

The inverled harmonic oscillator (IHO) is an unslable system, where Ihe Hamiltonian is given as 

1 


H = ^- 

2m 2 ^ 


(Al) 


In fhis system, Ihe classical Irajecfories are nol reslricled in a finile region, bul extends to infinilely 
large spalial and momenlum regions. While fhis unbounded nalure makes Ihe numerical calculation 
difficull, fhe analytic expression of Ihe HW entropy is known il5n . Then by comparing the numerical 
results with the analytic solution, we can examine the validity and the precision of the numerical 
methods. 


A.l. Analytic solution 

When the initial distribution of Wigner function is given by a Gaussian, 


fw{p,q;t = 0) = 2exp ( -^p 


U! 

h 


q 


(A2) 


the time evolution of the HW entropy is obtained analytically il5ll . Since the potential is quadratic, 
the semiclassical analysis is exact, and the time evolution of the Wigner function is calculated by 
solving the classical equation of motion. The HW entropy at time t is given as 11511 . 


Snwit) = log 

^/m 

2 

+1) 

(A3) 

A{t) = 2(cj/9cosh(2Af) -|-1 -|- 56') , 

(A4) 
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Fig. A1 Time dependence of the HW entropy of the inverted harmonic oscillator in the TP method. 


A.2. Numerical results with TP 

Figure lAll shows the time evolution of the HW entropy of IHO calculated in the TP method with 
Atp = 50 — 800. We find that the TP method can well descrihe the time evolution of the HW 
entropy at early times, and that numerical results show saturated behavior in later times. Since there 
exists a limiting value of S'hw in the TP method as discussed in Suhsec. l4^ we need to take a large 
number of App to describe a large amount of entropy production. It should be noted that numerical 
results converge in the limit of App —oo, and the converged result well describe the analytic result. 



t 


Fig. A2 Time dependence of the HW entropy of the inverted harmonic oscillator in the tsMC 
method. 


A. 3. Numerical results with tsMC 

Figure IA2I shows the time evolution of the HW entropy of IHO in the tsMC method with Ain = 
Nout = 100 and 1000. We find that numerical results are consistent with the analytic solution at 
early times t < 3, but that the numerical results tend to overestimate the analytic results and numer¬ 
ical errors become very large at later times. The large error would come from the poor overlap 
between the Wigner function and the coarse-graining Gaussian function at later time, which makes 
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importance sampling spoiled. On the other hand, the systematic overestimation may he due to the 
lack of sampling points in the Monte-Carlo integration in the logarithmic function (See Eq. (fT^ l. 
We note here that the statistical-error estimation is performed only for the Monte-Carlo integration 
outside the log. 
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